#ifndef RSIM_VECTOR_H
#define RSIM_VECTOR_H

/* --------------------------------------------------------------------------*
* File: Vector.h                                                             *
* Authors: Kishor Bhalerao                                                   *
* Email : kishor8dm@gmail.com                                                *
* Contributors:                                                              *
* Email:                                                                     *
* copyright (c) 2010 Authors.                                                *
*                                                                            *
* This program is free software: you can redistribute it and/or modify it    *
* under the terms of the GNU General Public License as published by the Free *
* Software Foundation, either version 3 of the License, or any later version.*
*                                                                            *
* This program is distributed in the hope that it will be useful, but        *
* WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
* for more details. (http://www.gnu.org/licenses/)                           *
* -------------------------------------------------------------------------- */

#include "CData.h"

#include <cstdlib>
#include <typeinfo>

using std::cout;
using std::endl;

namespace RSIM{

// template<class T>class Matrix_;

template <class T, int N, bool C=true> class Vector_;
template <class T, int N=0, bool C> class Vector_;
template <class T=double,int N, bool C> class Vector_;

template<class T, int N, bool C>
class Vector_:public cData<T,N>{
        protected:
                bool isColVec_;
        public:
                /** Create a vector of specified length */
                Vector_():cData<T,N>(),isColVec_(C){}
                
                /** Deep copy */
                template <class P,int Q, bool R>
                Vector_(const Vector_<P,Q,R>& V):cData<T,N>(),isColVec_(C){
                        #ifdef RSIM_SAFE
                        if(typeid(P)!=typeid(T)){
                                cout<<"Vector_(const Vector&): Attempt to convert Vector of one type into another\n";
                                exit(0);
                        }
                        this->resize(Q);
                        #endif
                        const P* ptr = V.Ptr();
                        for(int i=0;i<this->nele_;++i)
                                this->data_[i] = ptr[i];
                }
                                
                /** Deep copy */
//                 Vector_(const Matrix_<T>& M);
                
                /**
                Create a vector with all elements = val
                */
                Vector_(const T& val):cData<T,N>(val),isColVec_(C){}
                
                /** @return length of the vector */
                int dim()const{return this->nele_;}
                
                /** @return i^th element of the vector. 0 based indexing */
                const T& getElement(const int& i)const{return this->data_[i];}                                                                                
                
                /** toggles the isColVec variable. */
                void transpose(){isColVec_ = !isColVec_;}
                
                /** true if number of elements = 0 */
                bool isNull()const{return !(static_cast<bool>(this->nele_));}
                
                /** True if column vector */ 
                bool& isColVector(){return isColVec_;}
                
                /** True if column vector */
                bool isColVector()const{return isColVec_;}                                
                
                /** Returns a column vector*/
                Vector_<T> extractNElements(const int& n)const;
                
                void operator+=(const Vector_<T,N,C>& V);
                
                Vector_<T,N,C> operator-(const Vector_<T,N,C>& V)const;
                
                Vector_<T,N,C> operator+(const Vector_<T,N,C>& V)const;
                
//                 Vector_<T> operator-(const Matrix_<T>& M)const;
                
                /** vector x scalar */
                Vector_<T,N,C> operator*(const T& val)const;
                
                /** vector x matrix */
//                 Matrix_<T> operator*(const Matrix_<T>& M)const;
                
                /** vector x vector */
//                 Matrix_<T> operator*(const Vector_<T>& V)const;
                
                /** Transpose of a Matrix_ */
                Vector_<T,N,C> operator~(void) const;	
                
                /** @return writable reference */
                T& operator () (const int& i){return this->data_[i];}
                
                /** @return constant reference */
                T& operator () (const int& i) const{return this->data_[i];}
                
                /** @return writable reference */
                T& operator [] (const int& i){return this->data_[i];}
                
                /** @return constant reference */
                const T& operator [] (const int& i)const{return this->data_[i];}
                
                /** Deep copy */
                template<class P,int Q,bool R>
                void operator=(const Vector_<P,Q,R>& V){
                        #ifdef RSIM_SAFE
                        if(typeid(P)!=typeid(T)){
                                cout<<"Vector_(const Vector&): Attempt to convert Vector of one type into another\n";
                                exit(0);
                        }
                        this->resize(Q);
                        #endif                        
                        
                        const double *v = V.Ptr();
                        for(int i=0;i<this->nele_;++i){
                                this->data_[i] = v[i];
                        }
                }
                
                /** Deep copy */
//                 void operator=(const Matrix_<T>& M);
                
                template <class U,int Num,bool Co>
                friend std::ostream& operator<<(std::ostream& func, const Vector_<U,Num,Co>& V);
                
                template <class U,int Nf,bool Cf>
                friend Vector_<U,Nf,Cf> operator*(const double& val, const Vector_<U,Nf,Cf>& V);			
                
                template <class U,int Nf,bool Cf>
                friend std::istream& operator>>(std::istream& func, RSIM::Vector_<U>& V);

}; // class Vector_

////////////////////////////////////////////////////////////////////////

template<class T,int N,bool C>
Vector_<T> Vector_<T,N,C>::extractNElements(const int& n)const{
        Vector_<T> v;
        v.resize(n);
        double *ptr = v.wPtr();
        for(int i=0;i<n;++i)
                ptr[i]=this->data_[i];
        return v;
}

////////////////////////////////////////////////////////////////////////

template<class T,int N, bool C>
inline Vector_<T,N,C> Vector_<T,N,C>::operator~(void) const {
        Vector_<T,N,C> op(*this);
        op.isColVector() = !isColVec_;
        return op;
}

//////////////////////////////////////////////////////////////////////

template<class T, int N, bool C>
Vector_<T,N,C> Vector_<T,N,C>::operator+(const Vector_<T,N,C>& V)const{    
        Vector_<T,N,C> op;
        const T *v = V.Ptr();
        T *opvec = op.wPtr();
        for(int i=0;i<this->nele_;++i)
                opvec[i]=this->data_[i]+v[i];
        return op;
}

////////////////////////////////////////////////////////////////////////

template<class T, int N, bool C>
void Vector_<T,N,C>::operator+=(const Vector_<T,N,C>& V){                
        const T *v = V.Ptr();        
        for(int i=0;i<this->nele_;++i)
                this->data_[i] += v[i];        
}

////////////////////////////////////////////////////////////////////////

template<class T, int N, bool C>
Vector_<T,N,C> Vector_<T,N,C>::operator-(const Vector_<T,N,C>& V)const{        
        Vector_<T,N,C> op;
        const T *v = V.Ptr();
        T *opvec = op.wPtr();
        for(int i=0;i<this->nele_;++i)
                opvec[i]=this->data_[i]-v[i];
        return op;
}

////////////////////////////////////////////////////////////////////////
// 
// template<class T>
// Vector_<T> Vector_<T>::operator+(const Matrix_<T>& M)const{
//         int row = M.getNRow();
//         int col = M.getNCol();
// 
//         if(isColVec_ && (col==1)){
//                 Vector_<T> op(this->nele_);
//                 T *opvec = op.wPtr();
//                 
//                 const T *m = M.Ptr();
//                 for(int i=0;i<this->nele_;++i)
//                         opvec[i]=this->data_[i]+m[i];
//                 return op;
//         }
//         else if((!isColVec_) && (row == 1)){
//                 Vector_<T> op(this->nele_);
//                 op.isColVector() = false;
//                 T *opvec = op.wPtr();
//                 const T *m = M.Ptr();
//                 for(int i=0;i<this->nele_;++i)
//                         opvec[i] = this->data_[i]+m[i];
//                 return op;
//         }
//         cout<<"\nVector_ Vector_::operator+(const Matrix_<T>& M)const\n\
//         Inconsistent Dimensions in addition. exiting...\n";
//         exit(0);
// }
// 
// ////////////////////////////////////////////////////////////////////////
// 
// template<class T>
// Vector_<T> Vector_<T>::operator-(const Matrix_<T>& M)const{
//         int row = M.getNRow();
//         int col = M.getNCol();
//         if(isColVec_ && (col==1)){
//                 Vector_<T> op(this->nele_);
//                 T *opvec = op.wPtr();
//                 
//                 const T *m = M.Ptr();
//                 for(int i=0;i<this->nele_;++i)
//                         opvec[i]=this->data_[i]-m[i];
//                 return op;
//         }
//         else if(!isColVec_ && (row == 1)){
//                 Vector_<T> op(this->nele_);
//                 op.isColVector() = false;
//                 T *opvec = op.wPtr();
//                 const T *m = M.Ptr();
//                 for(int i=0;i<this->nele_;++i)
//                         opvec[i] = this->data_[i]-m[i];
//                 return op;
//         }
//         cout<<"\nVector_ Vector_::operator-(const Matrix_<T>& M)const\n\
//         Inconsistent Dimensions in addition. exiting...\n";
//         exit(0);
// }
// 
////////////////////////////////////////////////////////////////////////

template<class T,int N,bool C>
Vector_<T,N,C> Vector_<T,N,C>::operator*(const T& val)const{
        if(!this->nele_){
                std::cout<<"\nVector_ Vector_::operator*(const T& val)const\n\
                Attempt to use zero vector. About to exit...\n";
                exit(0);
        }
        Vector_<T,N,C> op(val);        
        
        T *opvec = op.wPtr();
        for(int i=0;i<N;i++){
                opvec[i]=this->data_[i]*val;
        }
        return op;
}

// ////////////////////////////////////////////////////////////////////////
// 
// template<class T>
// Matrix_<T> Vector_<T>::operator*(const Matrix_<T>& M)const{
//         int Arows, Acols, Brows, Bcols;
//         if (isColVec_){
//                 Arows = this->nele_;
//                 Acols = 1;
//         }
//         else{
//                 Arows = 1;
//                 Acols = this->nele_;
//         }
//         Brows = M.getNRow();
//         Bcols = M.getNCol();
//         #ifdef RSIM_SAFE
//         if(Acols!=Brows){
//                 cout<<"\nMatrix_<T> Vector_::operator*(const Matrix_<T>& M)const\n\
//                 Incompatible dimensions. About to exit...\n";
//                 exit(0);
//         }
//         #endif
//         Matrix_<T> op(Arows,Bcols);
//         T *opmat = op.wPtr();
//         const T *m = M.Ptr();
//         
//         T sum;
//         for(int i = 0;i<Arows;i++){
//                 for(int j=0;j<Bcols;j++){
//                         sum=0.0;
//                         for(int k=0;k<Acols;k++)
//                                 sum = sum+this->data_[i*Acols+k]*m[k*Bcols+j];
//                         opmat[i*Bcols+j]=sum;
//                 }
//         }
//         return op;
// }
// 
// ////////////////////////////////////////////////////////////////////////
// 
// template<class T>
// Matrix_<T> Vector_<T>::operator*(const Vector_<T>& V)const{
//         int Arows, Acols, Brows, Bcols;
//         if (isColVec_){
//                 Arows = this->nele_;
//                 Acols = 1;
//         }
//         else{
//                 Arows = 1;
//                 Acols = this->nele_;
//         }
//         if (V.isColVector()){
//                 Brows = V.getLength();
//                 Bcols = 1;
//         }
//         else{
//                 Brows = 1;
//                 Bcols = V.getLength();
// 
//         }	
//         #ifdef RSIM_SAFE
//         if(Acols!=Brows){
//                 cout<<"\nMatrix_<T> Vector_::operator*(const Vector_& V)const\n\
//                         Incompatible Matrix_<T> dimensions. About to exit...\n";
//                 exit(0);
//         }
//         #endif
//         Matrix_<T> op(Arows,Bcols);
//         T *opmat = op.wPtr();
//         const T *v2 = V.Ptr();
//         T sum;
//         for(int i = 0;i<Arows;i++){
//                 for(int j=0;j<Bcols;j++){
//                         sum=0.0;
//                         for(int k=0;k<Acols;k++)
//                                 sum = sum+this->data_[i+k]*v2[k+j];
//                         opmat[i*Bcols+j]=sum;
//                 }
//         }
//         return op;
// }
// 
// ////////////////////////////////////////////////////////////////////////
// 
// template<class T>
// void Vector_<T>::operator=(const Matrix_<T>& M){
//         int r = M.getNRow();
//         int c = M.getNCol();
//         if( !((r-1)*(c-1)) && !(this->nele_-r*c)){
//                 if(!(r-1))
//                         isColVec_ = false;
//                 else
//                         isColVec_ = true;
//                 const T *ptr = M.Ptr();
//                 for(int i=0;i<this->nele_;++i)
//                         this->data_[i] = ptr[i];
//         }
//         else if( !((r-1)*(c-1)) && !this->nele_){
//                 if(!(r-1))
//                         isColVec_ = false;
//                 else
//                         isColVec_ = true;
//                 this->setLength(r*c);
//                 const T *ptr = M.Ptr();
//                 for(int i=0;i<this->nele_;++i)
//                         this->data_[i] = ptr[i];
//         }
//         else{
//                 cout<<"\nvoid Vector_::operator=(const Matrix_<T>& M)\n\
//                 Incompatible dimensions, exiting... \n";
//                 exit(0);
//         }
// }
// 
// ////////////////////////////////////////////////////////////////////////
// 
// template<class T>
// Vector_<T>::Vector_(const Matrix_<T>& M):cData<T>(),isColVec_(true){
//         int r = M.getNRow();
//         int c = M.getNCol();
//         if( !((r-1)*(c-1)) ){
//                 if(!(r-1))
//                         isColVec_ = false;
//                 this->setLength(r*c);
//                 const T *ptr = M.Ptr();
//                 for(int i=0;i<this->nele_;++i)
//                         this->data_[i] = ptr[i];
//         }
//         else{
//                 cout<<"\nVector_::Vector_(const Matrix_<T>& M)\n\
//                 either row or col of the matrix M should be 1. Exiting... \n";
//                 exit(0);
//         }
// }
// ////////////////////////////////////////////////////////////////////////
// 
template<class T,int N, bool C>
std::ostream& operator<<(std::ostream& func, const Vector_<T,N,C>& V){ 
        if(V.isColVector())
                func<<"~[";
        else		
                func<<"[";
        for(int i=0;i<V.getLength();i++){
                if(i)
                        func<<", ";
                func<<V.getElement(i);
        }
        func<<"]";
        return func;
}
// 
////////////////////////////////////////////////////////////////////////

template <class T,int N,bool C>
Vector_<T,N,C> operator*(const double& val, const Vector_<T,N,C>& V){
        return (V*val);
}

////////////////////////////////////////////////////////////////////////

template <class T,int N,bool C>
std::istream& operator>>(std::istream& func, Vector_<T>& V){
        for(int i=0;i<V.getLength();i++){
                func>>V(i); 
        }
        return func; 
}

////////////////////////////////////////////////////////////////////////
}// namespace RSIM

#endif
